Analysis date: 2023-06-19

Depends on

DIPG_FirstBatch_DataProcessing Script

load("../Data/Cache/Xenografts_Batch1_2_DataProcessing.RData")

TODO

  • Do differential abudance analysis for prep batch and mass spec run

Setup

Load libraries and functions

Analysis

DEP

Tyrosine all

Each condition vs ctrl

Test_PhosphoData(pY_Set1_form, comparison = "E", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test
## Testing non-log transformed changes with t-test and adjusting p-value using FDR. Significance level = 0.05
data_diff_E_vs_ctrl_pY <- test_diff(pY_se_Set1, type="manual", test = "E_vs_ctrl")
## Tested contrasts: E_vs_ctrl
dep_E_vs_ctrl_pY <- add_rejections_SH(data_diff_E_vs_ctrl_pY, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_E_vs_ctrl_pY, contrast = "E_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pY") 
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Return_DEP_Hits_Plots(data = pY_Set1_form, dep_E_vs_ctrl_pY, comparison = "E_vs_ctrl_diff")
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(comparison)
## 
##   # Now:
##   data %>% select(all_of(comparison))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## 'select()' returned 1:1 mapping between keys and columns
## Loading required namespace: reactome.db
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                            pathway      pval
## 1: A tetrasaccharide linker sequence is required for GAG synthesis 0.3219316
## 2:                                       ABC transporter disorders 0.3501006
## 3:                          ABC-family proteins mediated transport 0.3501006
## 4:                       ADP signalling through P2Y purinoceptor 1 0.2273642
## 5:                                           ALK mutants bind TKIs 0.6953125
## 6:               APC/C-mediated degradation of cell cycle proteins 0.9683698
##         padj    log2err         ES        NES size leadingEdge
## 1: 0.8772981 0.10473282  0.8430233  1.1431478    1        6385
## 2: 0.8772981 0.09957912  0.8255814  1.1194965    1        5687
## 3: 0.8772981 0.09957912  0.8255814  1.1194965    1        5687
## 4: 0.8772981 0.12814292  0.8895349  1.2062180    1        1432
## 5: 0.9726046 0.06143641 -0.6569767 -0.8823673    1        1213
## 6: 0.9812774 0.05617666  0.3859649  0.6052109    2        5687

## Note: Row-scaling applied for this heatmap

Test_PhosphoData(pY_Set1_form, comparison = "EC", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test
## Testing non-log transformed changes with t-test and adjusting p-value using FDR. Significance level = 0.05
data_diff_EC_vs_ctrl_pY <- test_diff(pY_se_Set1, type="manual", test = "EC_vs_ctrl")
## Tested contrasts: EC_vs_ctrl
dep_EC_vs_ctrl_pY <- add_rejections_SH(data_diff_EC_vs_ctrl_pY, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EC_vs_ctrl_pY, contrast = "EC_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pY") 
Return_DEP_Hits_Plots(data = pY_Set1_form, dep_EC_vs_ctrl_pY, comparison = "EC_vs_ctrl_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                            pathway      pval
## 1: A tetrasaccharide linker sequence is required for GAG synthesis 0.6124031
## 2:                                       ABC transporter disorders 0.7926357
## 3:                          ABC-family proteins mediated transport 0.7926357
## 4:                       ADP signalling through P2Y purinoceptor 1 0.2984496
## 5:                                           ALK mutants bind TKIs 0.6404040
## 6:               APC/C-mediated degradation of cell cycle proteins 0.6906077
##       padj    log2err         ES        NES size leadingEdge
## 1: 0.97396 0.06720651  0.6918605  0.9340944    1        6385
## 2: 0.97396 0.05490737  0.5813953  0.7849533    1        5687
## 3: 0.97396 0.05490737  0.5813953  0.7849533    1        5687
## 4: 0.97396 0.10714024  0.8662791  1.1695804    1        1432
## 5: 0.97396 0.06705126 -0.6744186 -0.9142164    1        1213
## 6: 0.97396 0.05896945 -0.5518113 -0.8639961    2         983
## Warning in min(screen_pval05_neg[, logFcColStr]): no non-missing arguments to
## min; returning Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(cs1s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs1s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in min(cs2s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs2s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in min(cs3s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs3s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Note: Row-scaling applied for this heatmap

Plot_Enrichment_Single_Pathway(dep_EC_vs_ctrl_pY, comparison = "EC_vs_ctrl_diff", 
                               pw = "Epigenetic regulation of gene expression")
Test_PhosphoData(pY_Set1_form, comparison = "EBC", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test
## Testing non-log transformed changes with t-test and adjusting p-value using FDR. Significance level = 0.05
data_diff_EBC_vs_ctrl_pY <- test_diff(pY_se_Set1, type="manual", test = "EBC_vs_ctrl")
## Tested contrasts: EBC_vs_ctrl
dep_EBC_vs_ctrl_pY <- add_rejections_SH(data_diff_EBC_vs_ctrl_pY, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EBC_vs_ctrl_pY, contrast = "EBC_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pY")
Return_DEP_Hits_Plots(data = pY_Set1_form, dep_EBC_vs_ctrl_pY, comparison = "EBC_vs_ctrl_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                            pathway      pval
## 1: A tetrasaccharide linker sequence is required for GAG synthesis 0.4151329
## 2:                                       ABC transporter disorders 0.9213052
## 3:                          ABC-family proteins mediated transport 0.9213052
## 4:                       ADP signalling through P2Y purinoceptor 1 0.8875256
## 5:                                           ALK mutants bind TKIs 0.2053743
## 6:               APC/C-mediated degradation of cell cycle proteins 0.4480409
##         padj    log2err         ES        NES size leadingEdge
## 1: 0.9374656 0.09054289 -0.7790698 -1.0625762    1        6385
## 2: 0.9646724 0.04754342  0.5348837  0.7313084    1        5687
## 3: 0.9646724 0.04754342  0.5348837  0.7313084    1        5687
## 4: 0.9646724 0.05216303 -0.5581395 -0.7612486    1        1432
## 5: 0.8992707 0.13214726  0.8779070  1.2002996    1        1213
## 6: 0.9374656 0.07647671 -0.6679829 -1.0309284    2         983

## Note: Row-scaling applied for this heatmap

EC vs E

Test_PhosphoData(pY_Set1_form, comparison = "EC", comparison_base = "E") %>% GGPlotly_Volcano_Test
## Testing non-log transformed changes with t-test and adjusting p-value using FDR. Significance level = 0.05
data_diff_EC_vs_E_pY <- test_diff(pY_se_Set1, type = "manual", 
                              test = c("EC_vs_E"))
## Tested contrasts: EC_vs_E
dep_EC_vs_E_pY <- add_rejections_SH(data_diff_EC_vs_E_pY, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EC_vs_E_pY, contrast = "EC_vs_E",  add_names = TRUE, additional_title = "pY", proteins_of_interest = "EGFR")
Return_DEP_Hits_Plots(data = pY_Set1_form, dep_EC_vs_E_pY, comparison = "EC_vs_E_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns

##                                                            pathway      pval
## 1: A tetrasaccharide linker sequence is required for GAG synthesis 0.2920696
## 2:                                       ABC transporter disorders 0.1566731
## 3:                          ABC-family proteins mediated transport 0.1566731
## 4:                       ADP signalling through P2Y purinoceptor 1 0.8142857
## 5:                                           ALK mutants bind TKIs 0.1959184
## 6:               APC/C-mediated degradation of cell cycle proteins 0.3647059
##         padj    log2err         ES        NES size leadingEdge
## 1: 0.6828670 0.10839426 -0.8430233 -1.1437044    1        6385
## 2: 0.6329398 0.15419097 -0.9127907 -1.2383558    1        5687
## 3: 0.6329398 0.15419097 -0.9127907 -1.2383558    1        5687
## 4: 0.9108626 0.05605959  0.5930233  0.7873775    1        1432
## 5: 0.6486362 0.14040624  0.9069767  1.2042244    1        1213
## 6: 0.7323529 0.10672988 -0.7251462 -1.1327565    2    5687,983
#data_results <- get_df_long(dep)

EBC vs EC

Test_PhosphoData(pY_Set1_form, comparison = "EBC", comparison_base = "EC") %>% GGPlotly_Volcano_Test
## Testing non-log transformed changes with t-test and adjusting p-value using FDR. Significance level = 0.05
data_diff_EBC_vs_EC_pY <- test_diff(pY_se_Set1, type = "manual", 
                              test = c("EBC_vs_EC"))
## Tested contrasts: EBC_vs_EC
dep_EBC_vs_EC_pY <- add_rejections_SH(data_diff_EBC_vs_EC_pY, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EBC_vs_EC_pY, contrast = "EBC_vs_EC",  add_names = TRUE, additional_title = "pY")
Return_DEP_Hits_Plots(data = pY_Set1_form, dep_EBC_vs_EC_pY, comparison = "EBC_vs_EC_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns

##                                                            pathway      pval
## 1: A tetrasaccharide linker sequence is required for GAG synthesis 0.1011673
## 2:                                       ABC transporter disorders 0.9533074
## 3:                          ABC-family proteins mediated transport 0.9533074
## 4:                       ADP signalling through P2Y purinoceptor 1 0.1595331
## 5:                                           ALK mutants bind TKIs 0.3385214
## 6:               APC/C-mediated degradation of cell cycle proteins 0.8342541
##         padj    log2err         ES        NES size leadingEdge
## 1: 0.4419371 0.19578900 -0.9476744 -1.2667050    1        6385
## 2: 0.9748400 0.04660151 -0.5232558 -0.6994077    1        5687
## 3: 0.9748400 0.04660151 -0.5232558 -0.6994077    1        5687
## 4: 0.5884998 0.15315881 -0.9186047 -1.2278491    1        1432
## 5: 0.8265848 0.09957912 -0.8313953 -1.1112811    1        1213
## 6: 0.9748400 0.05019343 -0.5263158 -0.7832006    2    983,5687
#data_results <- get_df_long(dep)

Serine/Threonine fractionated

Each condition vs ctrl

#Test_PhosphoData(pST_Set1_form, comparison = "E", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test

data_diff_E_vs_ctrl_pST <- test_diff(pST_se_Set1, type="manual", test = "E_vs_ctrl")
## Tested contrasts: E_vs_ctrl
dep_E_vs_ctrl_pST <- add_rejections_SH(data_diff_E_vs_ctrl_pST, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_E_vs_ctrl_pST, contrast = "E_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pST") 
Return_DEP_Hits_Plots(data = pST_Set1_form, dep_E_vs_ctrl_pST, comparison = "E_vs_ctrl_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                    pathway      pval      padj
## 1:                                  2-LTR circle formation 0.4884211 0.8732408
## 2:                               ABC transporter disorders 0.4390681 0.8732408
## 3:                  ABC-family proteins mediated transport 0.9052632 0.9743672
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.3179298 0.8732408
## 5:               ADP signalling through P2Y purinoceptor 1 0.1797323 0.7981031
## 6:               AKT phosphorylates targets in the cytosol 0.3315789 0.8732408
##       log2err         ES        NES size             leadingEdge
## 1: 0.08312913 -0.5815003 -1.0045958    3              11168,2547
## 2: 0.08020234  0.5617310  1.0556435    5     5684,5707,9491,5708
## 3: 0.04424074  0.3164369  0.6180787    6     5684,5707,4363,9491
## 4: 0.10027911  0.6518453  1.1434323    4               5576,5573
## 5: 0.14205664 -0.9150889 -1.2188049    1                    5321
## 6: 0.09466462  0.5827565  1.1382658    6 7249,84335,207,572,2931

## Note: Row-scaling applied for this heatmap

#Test_PhosphoData(pST_Set1_form, comparison = "EC", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test

data_diff_EC_vs_ctrl_pST <- test_diff(pST_se_Set1, type="manual", test = "EC_vs_ctrl")
## Tested contrasts: EC_vs_ctrl
dep_EC_vs_ctrl_pST <- add_rejections_SH(data_diff_EC_vs_ctrl_pST, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EC_vs_ctrl_pST, contrast = "EC_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pST") 
Return_DEP_Hits_Plots(data = pST_Set1_form, dep_EC_vs_ctrl_pST, comparison = "EC_vs_ctrl_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                    pathway      pval      padj
## 1:                                  2-LTR circle formation 0.7790262 0.9423051
## 2:                               ABC transporter disorders 0.3752711 0.8362723
## 3:                  ABC-family proteins mediated transport 0.6033403 0.8723999
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.2838710 0.7967630
## 5:               ADP signalling through P2Y purinoceptor 1 0.3073930 0.7967630
## 6:               AKT phosphorylates targets in the cytosol 0.6450939 0.8848830
##       log2err         ES        NES size               leadingEdge
## 1: 0.05412006 -0.4836301 -0.8035502    3           11168,2547,3159
## 2: 0.09992770  0.5650003  1.0950764    5 5707,9491,5708,55788,5684
## 3: 0.07165274  0.4391696  0.8874789    6            5707,9491,4363
## 4: 0.11724972  0.6561695  1.2011347    4                 5576,5573
## 5: 0.10552094  0.8531268  1.1330559    1                      5321
## 6: 0.06831109  0.4199677  0.8486755    6           7249,84335,2931

## Note: Row-scaling applied for this heatmap

Plot_Enrichment_Single_Pathway(dep_EC_vs_ctrl_pST, comparison = "EC_vs_ctrl_diff", 
                               pw = "Epigenetic regulation of gene expression")
#Test_PhosphoData(pST_Set1_form, comparison = "EBC", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test

data_diff_EBC_vs_ctrl_pST <- test_diff(pST_se_Set1, type="manual", test = "EBC_vs_ctrl")
## Tested contrasts: EBC_vs_ctrl
dep_EBC_vs_ctrl_pST <- add_rejections_SH(data_diff_EBC_vs_ctrl_pST, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EBC_vs_ctrl_pST, contrast = "EBC_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pST")
Return_DEP_Hits_Plots(data = pST_Set1_form, dep_EBC_vs_ctrl_pST, comparison = "EBC_vs_ctrl_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns

##                                                    pathway      pval      padj
## 1:                                  2-LTR circle formation 0.1364366 0.7522345
## 2:                               ABC transporter disorders 0.1638889 0.7539830
## 3:                  ABC-family proteins mediated transport 0.8616541 0.9691266
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.6256831 0.8890290
## 5:               ADP signalling through P2Y purinoceptor 1 0.2534930 0.7664718
## 6:               AKT phosphorylates targets in the cytosol 0.8160237 0.9509037
##       log2err         ES        NES size               leadingEdge
## 1: 0.15016980 -0.7587593 -1.3124585    3           2547,3159,11168
## 2: 0.18302394  0.5767683  1.3037564    5 9491,5707,55788,5708,5684
## 3: 0.03943665 -0.3321983 -0.6955622    6                        23
## 4: 0.08383611  0.4209305  0.8668766    4                 5576,5573
## 5: 0.11988785 -0.8703385 -1.1636827    1                      5321
## 6: 0.07417590  0.3029458  0.7300750    6            84335,207,1026

EC vs E

#Test_PhosphoData(pST_Set1_form, comparison = "EC", comparison_base = "E") %>% GGPlotly_Volcano_Test

data_diff_EC_vs_E_pST <- test_diff(pST_se_Set1, type = "manual", 
                              test = c("EC_vs_E"))
## Tested contrasts: EC_vs_E
dep_EC_vs_E_pST <- add_rejections_SH(data_diff_EC_vs_E_pST, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EC_vs_E_pST, contrast = "EC_vs_E",  add_names = TRUE, additional_title = "pST", proteins_of_interest = "EGFR")
Return_DEP_Hits_Plots(data = pST_Set1_form, dep_EC_vs_E_pST, comparison = "EC_vs_E_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                    pathway        pval
## 1:                                  2-LTR circle formation 0.427792916
## 2:                               ABC transporter disorders 0.493525180
## 3:                  ABC-family proteins mediated transport 0.779944290
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.503681885
## 5:               ADP signalling through P2Y purinoceptor 1 0.013641932
## 6:               AKT phosphorylates targets in the cytosol 0.008241885
##         padj    log2err         ES        NES size        leadingEdge
## 1: 0.9885312 0.10592029  0.5614847  1.0256258    3         11168,2547
## 2: 0.9885312 0.06321912 -0.5018511 -1.0017106    5               5684
## 3: 0.9885312 0.04049348 -0.3692819 -0.7751279    6               5684
## 4: 0.9885312 0.06335970 -0.5293804 -0.9870411    4     5576,5577,5573
## 5: 0.6825367 0.38073040  0.9913941  1.3389465    1               5321
## 6: 0.6479679 0.38073040 -0.7830913 -1.6437198    6 572,84335,7249,207

## Note: Row-scaling applied for this heatmap

#data_results <- get_df_long(dep)

EBC vs EC

#Test_PhosphoData(pST_Set1_form, comparison = "EBC", comparison_base = "EC") %>% GGPlotly_Volcano_Test

data_diff_EBC_vs_EC_pST <- test_diff(pST_se_Set1, type = "manual", 
                              test = c("EBC_vs_EC"))
## Tested contrasts: EBC_vs_EC
dep_EBC_vs_EC_pST <- add_rejections_SH(data_diff_EBC_vs_EC_pST, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EBC_vs_EC_pST, contrast = "EBC_vs_EC",  add_names = TRUE, additional_title = "pST")
Return_DEP_Hits_Plots(data = pST_Set1_form, dep_EBC_vs_EC_pST, comparison = "EBC_vs_EC_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                    pathway       pval      padj
## 1:                                  2-LTR circle formation 0.14456036 0.7828541
## 2:                               ABC transporter disorders 0.60357143 0.8944963
## 3:                  ABC-family proteins mediated transport 0.36122178 0.8210576
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.69121813 0.9042184
## 5:               ADP signalling through P2Y purinoceptor 1 0.05642023 0.6434396
## 6:               AKT phosphorylates targets in the cytosol 0.30923695 0.8202010
##       log2err         ES        NES size          leadingEdge
## 1: 0.13959967 -0.7489948 -1.2676540    3      3159,11168,2547
## 2: 0.10135074  0.3812536  0.8866448    5 55788,5708,9491,5684
## 3: 0.07473852 -0.5384335 -1.0991481    6              23,4363
## 4: 0.04678830 -0.4540862 -0.8372966    4            5576,5577
## 5: 0.26635066 -0.9787722 -1.2999630    1                 5321
## 6: 0.15851411  0.4545464  1.1251306    6   1026,572,207,84335

## Note: Row-scaling applied for this heatmap

#data_results <- get_df_long(dep)

Serine/Threonine one shot

Each condition vs ctrl

#Test_PhosphoData(pST_Set1_oneshot_form, comparison = "E", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test

data_diff_E_vs_ctrl_pST_oneshot <- test_diff(pST_se_Set1_oneshot, type="manual", test = "E_vs_ctrl")
## Tested contrasts: E_vs_ctrl
dep_E_vs_ctrl_pST_oneshot <- add_rejections_SH(data_diff_E_vs_ctrl_pST_oneshot, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_E_vs_ctrl_pST_oneshot, contrast = "E_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pST") 
Return_DEP_Hits_Plots(data = pST_Set1_oneshot_form, dep_E_vs_ctrl_pST_oneshot, comparison = "E_vs_ctrl_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                    pathway       pval      padj
## 1:                                  2-LTR circle formation 0.55458515 0.8677316
## 2:                               ABC transporter disorders 0.29690722 0.8188248
## 3:                  ABC-family proteins mediated transport 0.69126214 0.9109398
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.07572816 0.6744680
## 5:               AKT phosphorylates targets in the cytosol 0.82785300 0.9580505
## 6:                                   ALK mutants bind TKIs 0.17361111 0.7510412
##       log2err         ES        NES size          leadingEdge
## 1: 0.07829552 -0.5072386 -0.9032856    3           11168,3159
## 2: 0.11146267  0.8545455  1.1470004    1                 5684
## 3: 0.06143641  0.5572177  0.8477599    2              5684,23
## 4: 0.22798720  0.9016990  1.3718593    2            5576,5573
## 5: 0.05280500 -0.5915152 -0.7872001    1                 7249
## 6: 0.13725078  0.6578377  1.2853796    5 5573,27436,1213,6801

## Note: Row-scaling applied for this heatmap

#Test_PhosphoData(pST_Set1_oneshot_form, comparison = "EC", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test

data_diff_EC_vs_ctrl_pST_oneshot <- test_diff(pST_se_Set1_oneshot, type="manual", test = "EC_vs_ctrl")
## Tested contrasts: EC_vs_ctrl
dep_EC_vs_ctrl_pST_oneshot <- add_rejections_SH(data_diff_EC_vs_ctrl_pST_oneshot, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EC_vs_ctrl_pST_oneshot, contrast = "EC_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pST") 
Return_DEP_Hits_Plots(data = pST_Set1_oneshot_form, dep_EC_vs_ctrl_pST_oneshot, comparison = "EC_vs_ctrl_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                    pathway       pval      padj
## 1:                                  2-LTR circle formation 0.69052632 0.9827737
## 2:                               ABC transporter disorders 0.69488189 0.9827737
## 3:                  ABC-family proteins mediated transport 0.42959002 0.9157638
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.13151927 0.7502400
## 5:               AKT phosphorylates targets in the cytosol 0.95669291 0.9919874
## 6:                                   ALK mutants bind TKIs 0.09368192 0.7099567
##       log2err         ES        NES size               leadingEdge
## 1: 0.06538342  0.4572435  0.8063100    3           3159,11168,7518
## 2: 0.06184060 -0.6545455 -0.8774269    1                      5684
## 3: 0.08108021 -0.6553398 -1.0345104    2                   23,5684
## 4: 0.18470647  0.8662491  1.3300696    2                 5576,5573
## 5: 0.04697587 -0.5175758 -0.6938172    1                      7249
## 6: 0.21654284  0.7422244  1.5241204    5 1213,27436,6801,5573,4869

## Note: Row-scaling applied for this heatmap

Plot_Enrichment_Single_Pathway(dep_EC_vs_ctrl_pST_oneshot, comparison = "EC_vs_ctrl_diff", 
                               pw = "Epigenetic regulation of gene expression")
#Test_PhosphoData(pST_Set1_oneshot_form, comparison = "EBC", comparison_base = "ctrl") %>% GGPlotly_Volcano_Test

data_diff_EBC_vs_ctrl_pST_oneshot <- test_diff(pST_se_Set1_oneshot, type="manual", test = "EBC_vs_ctrl")
## Tested contrasts: EBC_vs_ctrl
dep_EBC_vs_ctrl_pST_oneshot <- add_rejections_SH(data_diff_EBC_vs_ctrl_pST_oneshot, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EBC_vs_ctrl_pST_oneshot, contrast = "EBC_vs_ctrl", 
                 add_names = TRUE,
                additional_title = "pST")
Return_DEP_Hits_Plots(data = pST_Set1_oneshot_form, dep_EBC_vs_ctrl_pST_oneshot, comparison = "EBC_vs_ctrl_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                    pathway        pval
## 1:                                  2-LTR circle formation 0.826606876
## 2:                               ABC transporter disorders 0.773195876
## 3:                  ABC-family proteins mediated transport 0.533980583
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.006799441
## 5:               AKT phosphorylates targets in the cytosol 0.481624758
## 6:                                   ALK mutants bind TKIs 0.020092383
##         padj    log2err         ES        NES size               leadingEdge
## 1: 0.9402403 0.04107133 -0.4176059 -0.7462942    3                3159,11168
## 2: 0.9194048 0.05896945 -0.6060606 -0.8056837    1                      5684
## 3: 0.9194048 0.06508776 -0.6067961 -0.9614707    2                   23,5684
## 4: 0.2871996 0.40701792  0.9538835  1.5451380    2                 5576,5573
## 5: 0.9194048 0.07934350  0.7709091  1.0200925    1                      7249
## 6: 0.3168414 0.35248786  0.7259440  1.6927374    5 5573,27436,6801,4869,1213
## Warning in max(screen_pval05_pos[, logFcColStr]): no non-missing arguments to
## max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(cs1s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs1s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in min(cs2s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs2s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in min(cs3s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs3s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

EC vs E

#Test_PhosphoData(pST_Set1_oneshot_form, comparison = "EC", comparison_base = "E") %>% GGPlotly_Volcano_Test

data_diff_EC_vs_E_pST_oneshot <- test_diff(pST_se_Set1_oneshot, type = "manual", 
                              test = c("EC_vs_E"))
## Tested contrasts: EC_vs_E
dep_EC_vs_E_pST_oneshot <- add_rejections_SH(data_diff_EC_vs_E_pST_oneshot, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EC_vs_E_pST_oneshot, contrast = "EC_vs_E",  add_names = TRUE, additional_title = "pST", proteins_of_interest = "EGFR")
Return_DEP_Hits_Plots(data = pST_Set1_oneshot_form, dep_EC_vs_E_pST_oneshot, comparison = "EC_vs_E_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
## Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
## minSize, : For some of the pathways the P-values were likely overestimated. For
## such pathways log2err is set to NA.
##                                                    pathway       pval      padj
## 1:                                  2-LTR circle formation 0.30573248 0.9954594
## 2:                               ABC transporter disorders 0.05182342 0.6249294
## 3:                  ABC-family proteins mediated transport 0.08794788 0.8195144
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.92182410 0.9954594
## 5:               AKT phosphorylates targets in the cytosol 0.95010395 0.9954594
## 6:                                   ALK mutants bind TKIs 0.95786517 0.9954594
##       log2err         ES        NES size leadingEdge
## 1: 0.14040624  0.6081637  1.1278061    3  3159,11168
## 2: 0.27650060 -0.9757576 -1.2980569    1        5684
## 3: 0.19189224 -0.8714427 -1.3446052    2        5684
## 4: 0.04000315 -0.4524340 -0.6980896    2        5573
## 5: 0.04979032  0.5200000  0.6932750    1        7249
## 6: 0.03120530 -0.3004947 -0.5875252    5  5573,27436
## Warning in max(screen_pval05_pos[, logFcColStr]): no non-missing arguments to
## max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(cs1s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs1s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in min(cs2s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs2s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in min(cs3s, na.rm = TRUE): no non-missing arguments to min; returning
## Inf
## Warning in max(cs3s, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Note: Row-scaling applied for this heatmap

#data_results <- get_df_long(dep)

EBC vs EC

#Test_PhosphoData(pST_Set1_oneshot_form, comparison = "EBC", comparison_base = "EC") %>% GGPlotly_Volcano_Test

data_diff_EBC_vs_EC_pST_oneshot <- test_diff(pST_se_Set1_oneshot, type = "manual", 
                              test = c("EBC_vs_EC"))
## Tested contrasts: EBC_vs_EC
dep_EBC_vs_EC_pST_oneshot <- add_rejections_SH(data_diff_EBC_vs_EC_pST_oneshot, alpha = 0.05, lfc = log2(1.2))
GGPlotly_Volcano(dep_EBC_vs_EC_pST_oneshot, contrast = "EBC_vs_EC",  add_names = TRUE, additional_title = "pST")
Return_DEP_Hits_Plots(data = pST_Set1_oneshot_form, dep_EBC_vs_EC_pST_oneshot, comparison = "EBC_vs_EC_diff")
## 'select()' returned 1:1 mapping between keys and columns
## 'select()' returned 1:many mapping between keys and columns
## 'select()' returned 1:1 mapping between keys and columns
##                                                    pathway      pval      padj
## 1:                                  2-LTR circle formation 0.1170370 0.5876732
## 2:                               ABC transporter disorders 0.8133874 0.9592315
## 3:                  ABC-family proteins mediated transport 0.7240829 0.9592315
## 4: ADORA2B mediated anti-inflammatory cytokines production 0.6794258 0.9592315
## 5:               AKT phosphorylates targets in the cytosol 0.2454361 0.7566075
## 6:                                   ALK mutants bind TKIs 0.6816380 0.9592315
##       log2err         ES        NES size leadingEdge
## 1: 0.15631240 -0.7530280 -1.3088694    3        3159
## 2: 0.05582647  0.5866667  0.7898633    1        5684
## 3: 0.05009229 -0.5267946 -0.8318367    2          23
## 4: 0.05302125 -0.5463292 -0.8626828    2        5576
## 5: 0.12325723  0.8824242  1.1880589    1        7249
## 6: 0.04424074 -0.4274849 -0.8502453    5        1213

## Note: Row-scaling applied for this heatmap

#data_results <- get_df_long(dep)

Motif enrichment analysis

7AA Sequences pST fractionated

all_pST_sites <- (rowData(dep_E_vs_ctrl_pST) %>% as_tibble()) %>% select( Annotated_Sequence, HGNC_Symbol, Master.Protein.Accessions) %>% unique

table(all_pST_sites$Master.Protein.Accessions %in% Phosphorylation_site_dataset$ACC_ID )
## 
## FALSE  TRUE 
##    15  3923
Phosphorylation_site_dataset_pST <- Phosphorylation_site_dataset %>% filter(ACC_ID %in% all_pST_sites$Master.Protein.Accessions)

# TODO: s/t in first place might have TMT and or phos
# TODO did not account for double phosphorylated peptides
extracted_trimmed_peptide <- function(peptide){
  str_extract(peptide, "s|t")
  
  if( nchar( sub("[ts].*", "", peptide )) >= 7 & 
      nchar( sub(".*[ts]", "", peptide )) >= 7){
    pep_trimmed <- str_extract(peptide, ".......[ts].......")
  }else if(nchar( sub("[ts].*", "", peptide )) < 7 &
           nchar( sub(".*[ts]", "", peptide )) < 7){
    pep_trimmed <- str_extract(peptide, ".*[ts].*")
  }else if(nchar( sub("[ts].*", "", peptide )) < 7){
    pep_trimmed <- str_extract(peptide, ".*[ts].......")
  }else if(nchar( sub(".*[ts]", "", peptide )) < 7){
    pep_trimmed <- str_extract(peptide, ".......[ts].*")
  }else{stop(paste("Could not extract trimmed sequenxce from", peptide )) }
  return(pep_trimmed)
}

all_pST_sites$trimmed_sequence <- sapply(all_pST_sites$Annotated_Sequence, extracted_trimmed_peptide, simplify = T)
 
get_plusminus_7AA_sequence <- function(peptide, protein){
  mod_peptide <- paste0(str_to_upper(sub("[ts].*", "", peptide) ),
         str_extract(peptide, "[ts]"),
         str_to_upper(sub(".*[ts]", "", peptide) ) )
  sequ <- Phosphorylation_site_dataset_pST %>% filter(ACC_ID == protein) %>%
    mutate( SEQUENCE = paste0(str_to_upper(str_sub(`SITE_+/-7_AA`, 1, 7)),
                               str_sub(`SITE_+/-7_AA`, 8,8),
                              paste0(str_to_upper(str_sub(`SITE_+/-7_AA`, 9, 15)))) ) %>% 
    filter(str_detect(SEQUENCE, mod_peptide)) %>% .$SEQUENCE
  if(length(sequ) == 0){
    return(NA)}else{return(sequ)}
}

sequences7AA <- sapply(1:nrow(all_pST_sites), function(i){
  sequ <- get_plusminus_7AA_sequence(peptide = all_pST_sites$trimmed_sequence[i],
                             protein = all_pST_sites$Master.Protein.Accessions[i]) %>% unique
  if(length(sequ) != 1){
    return(NA)}else{return(sequ)}
  }, simplify = T) 

table(is.na(sequences7AA) )
## 
## FALSE  TRUE 
##  3493   445
all_pST_sites$Sequence_7_AA <- sequences7AA

7AA Sequences pST one shot

all_pST_oneshot_sites <- (rowData(dep_E_vs_ctrl_pST_oneshot) %>% as_tibble()) %>% select( Annotated_Sequence, HGNC_Symbol, Master.Protein.Accessions) %>% unique

table(all_pST_oneshot_sites$Master.Protein.Accessions %in% Phosphorylation_site_dataset$ACC_ID )
## 
## FALSE  TRUE 
##     1  1446
Phosphorylation_site_dataset_pST_oneshot <- Phosphorylation_site_dataset %>% filter(ACC_ID %in% all_pST_oneshot_sites$Master.Protein.Accessions)

# TODO: s/t in first place might have TMT and or phos
# TODO did not account for double phosphorylated peptides
extracted_trimmed_peptide <- function(peptide){
  str_extract(peptide, "s|t")
  
  if( nchar( sub("[ts].*", "", peptide )) >= 7 & 
      nchar( sub(".*[ts]", "", peptide )) >= 7){
    pep_trimmed <- str_extract(peptide, ".......[ts].......")
  }else if(nchar( sub("[ts].*", "", peptide )) < 7 &
           nchar( sub(".*[ts]", "", peptide )) < 7){
    pep_trimmed <- str_extract(peptide, ".*[ts].*")
  }else if(nchar( sub("[ts].*", "", peptide )) < 7){
    pep_trimmed <- str_extract(peptide, ".*[ts].......")
  }else if(nchar( sub(".*[ts]", "", peptide )) < 7){
    pep_trimmed <- str_extract(peptide, ".......[ts].*")
  }else{stop(paste("Could not extract trimmed sequenxce from", peptide )) }
  return(pep_trimmed)
}

all_pST_oneshot_sites$trimmed_sequence <- sapply(all_pST_oneshot_sites$Annotated_Sequence, extracted_trimmed_peptide, simplify = T)
 
get_plusminus_7AA_sequence <- function(peptide, protein){
  mod_peptide <- paste0(str_to_upper(sub("[ts].*", "", peptide) ),
         str_extract(peptide, "[ts]"),
         str_to_upper(sub(".*[ts]", "", peptide) ) )
  sequ <- Phosphorylation_site_dataset_pST_oneshot %>% filter(ACC_ID == protein) %>%
    mutate( SEQUENCE = paste0(str_to_upper(str_sub(`SITE_+/-7_AA`, 1, 7)),
                               str_sub(`SITE_+/-7_AA`, 8,8),
                              paste0(str_to_upper(str_sub(`SITE_+/-7_AA`, 9, 15)))) ) %>% 
    filter(str_detect(SEQUENCE, mod_peptide)) %>% .$SEQUENCE
  if(length(sequ) == 0){
    return(NA)}else{return(sequ)}
}

sequences7AA <- sapply(1:nrow(all_pST_oneshot_sites), function(i){
  sequ <- get_plusminus_7AA_sequence(peptide = all_pST_oneshot_sites$trimmed_sequence[i],
                             protein = all_pST_oneshot_sites$Master.Protein.Accessions[i]) %>% unique
  if(length(sequ) != 1){
    return(NA)}else{return(sequ)}
  }, simplify = T) 

table(is.na(sequences7AA) )
## 
## FALSE  TRUE 
##  1295   152
all_pST_oneshot_sites$Sequence_7_AA <- sequences7AA

E vs ctrl

E_vs_ctrl_pST_7AA <- left_join( (rowData(dep_E_vs_ctrl_pST) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = E_vs_ctrl_diff,  
           p = E_vs_ctrl_p.adj) 
    %>% unique ),
  (all_pST_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


E_vs_ctrl_pST_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/E_vs_ctrl_pST.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

E vs ctrl on shot

E_vs_ctrl_pST_oneshot_7AA <- left_join( (rowData(dep_E_vs_ctrl_pST_oneshot) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = E_vs_ctrl_diff,  
           p = E_vs_ctrl_p.adj) 
    %>% unique ),
  (all_pST_oneshot_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


E_vs_ctrl_pST_oneshot_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/E_vs_ctrl_pST_oneshot.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

EC vs ctrl

EC_vs_ctrl_pST_7AA <- left_join( (rowData(dep_EC_vs_ctrl_pST) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = EC_vs_ctrl_diff,  
           p = EC_vs_ctrl_p.adj) 
    %>% unique ),
  (all_pST_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


EC_vs_ctrl_pST_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/EC_vs_ctrl_pST.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

EC vs ctrl on shot

EC_vs_ctrl_pST_oneshot_7AA <- left_join( (rowData(dep_EC_vs_ctrl_pST_oneshot) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = EC_vs_ctrl_diff,  
           p = EC_vs_ctrl_p.adj) 
    %>% unique ),
  (all_pST_oneshot_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


EC_vs_ctrl_pST_oneshot_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/EC_vs_ctrl_pST_oneshot.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

EC vs E

EC_vs_E_pST_7AA <- left_join( (rowData(dep_EC_vs_E_pST) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = EC_vs_E_diff,  
           p = EC_vs_E_p.adj) 
    %>% unique ),
  (all_pST_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


EC_vs_E_pST_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/EC_vs_E_pST.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

EC vs E on shot

EC_vs_E_pST_oneshot_7AA <- left_join( (rowData(dep_EC_vs_E_pST_oneshot) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = EC_vs_E_diff,  
           p = EC_vs_E_p.adj) 
    %>% unique ),
  (all_pST_oneshot_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


EC_vs_E_pST_oneshot_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/EC_vs_E_pST_oneshot.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

EBC vs ctrl

EBC_vs_ctrl_pST_7AA <- left_join( (rowData(dep_EBC_vs_ctrl_pST) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = EBC_vs_ctrl_diff,  
           p = EBC_vs_ctrl_p.adj) 
    %>% unique ),
  (all_pST_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


EBC_vs_ctrl_pST_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/EBC_vs_ctrl_pST.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

EBC vs ctrl on shot

EBC_vs_ctrl_pST_oneshot_7AA <- left_join( (rowData(dep_EBC_vs_ctrl_pST_oneshot) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = EBC_vs_ctrl_diff,  
           p = EBC_vs_ctrl_p.adj) 
    %>% unique ),
  (all_pST_oneshot_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


EBC_vs_ctrl_pST_oneshot_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/EBC_vs_ctrl_pST_oneshot.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

EBC vs EC

EBC_vs_EC_pST_7AA <- left_join( (rowData(dep_EBC_vs_EC_pST) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = EBC_vs_EC_diff,  
           p = EBC_vs_EC_p.adj) 
    %>% unique ),
  (all_pST_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


EBC_vs_EC_pST_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/EBC_vs_EC_pST.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

EBC vs EC on shot

EBC_vs_EC_pST_oneshot_7AA <- left_join( (rowData(dep_EBC_vs_EC_pST_oneshot) %>% as_tibble() %>% 
    select(annotation = ID, Annotated_Sequence, HGNC_Symbol, 
           fc = EBC_vs_EC_diff,  
           p = EBC_vs_EC_p.adj) 
    %>% unique ),
  (all_pST_oneshot_sites %>% select(Annotated_Sequence, Sequence_7_AA, HGNC_Symbol)),
  by=c("Annotated_Sequence", "HGNC_Symbol") ) %>% filter(!is.na(Sequence_7_AA)) %>% 
  mutate( peptide = str_to_upper(Sequence_7_AA) ) %>%
  select(annotation, peptide, fc, p) %>% as.data.frame()


EBC_vs_EC_pST_oneshot_7AA %>% select(peptide, fc, p)  %>% write.table(file = "For_motif_analysis/EBC_vs_EC_pST_oneshot.txt", quote = FALSE, row.names = F, col.names = F, sep = "\t")

Session Info

sessionInfo()
## R version 4.2.3 (2023-03-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] lubridate_1.9.2             forcats_1.0.0              
##  [3] stringr_1.5.0               dplyr_1.1.2                
##  [5] purrr_1.0.1                 readr_2.1.4                
##  [7] tidyr_1.3.0                 tibble_3.2.1               
##  [9] ggplot2_3.4.2               tidyverse_2.0.0            
## [11] mdatools_0.14.0             SummarizedExperiment_1.28.0
## [13] GenomicRanges_1.50.2        GenomeInfoDb_1.34.9        
## [15] MatrixGenerics_1.10.0       matrixStats_1.0.0          
## [17] DEP_1.20.0                  org.Hs.eg.db_3.16.0        
## [19] AnnotationDbi_1.60.2        IRanges_2.32.0             
## [21] S4Vectors_0.36.2            Biobase_2.58.0             
## [23] BiocGenerics_0.44.0         fgsea_1.24.0               
## 
## loaded via a namespace (and not attached):
##   [1] circlize_0.4.15        fastmatch_1.1-3        plyr_1.8.8            
##   [4] igraph_1.5.0           gmm_1.8                lazyeval_0.2.2        
##   [7] shinydashboard_0.7.2   crosstalk_1.2.0        BiocParallel_1.32.6   
##  [10] digest_0.6.31          foreach_1.5.2          htmltools_0.5.5       
##  [13] fansi_1.0.4            magrittr_2.0.3         memoise_2.0.1         
##  [16] cluster_2.1.4          doParallel_1.0.17      tzdb_0.4.0            
##  [19] limma_3.54.2           ComplexHeatmap_2.14.0  Biostrings_2.66.0     
##  [22] imputeLCMD_2.1         sandwich_3.0-2         timechange_0.2.0      
##  [25] colorspace_2.1-0       blob_1.2.4             xfun_0.39             
##  [28] crayon_1.5.2           RCurl_1.98-1.12        jsonlite_1.8.5        
##  [31] impute_1.72.3          zoo_1.8-12             iterators_1.0.14      
##  [34] glue_1.6.2             hash_2.2.6.2           gtable_0.3.3          
##  [37] zlibbioc_1.44.0        XVector_0.38.0         GetoptLong_1.0.5      
##  [40] DelayedArray_0.24.0    shape_1.4.6            scales_1.2.1          
##  [43] pheatmap_1.0.12        vsn_3.66.0             mvtnorm_1.2-2         
##  [46] DBI_1.1.3              Rcpp_1.0.10            plotrix_3.8-2         
##  [49] mzR_2.32.0             viridisLite_0.4.2      xtable_1.8-4          
##  [52] clue_0.3-64            reactome.db_1.82.0     bit_4.0.5             
##  [55] preprocessCore_1.60.2  sqldf_0.4-11           MsCoreUtils_1.10.0    
##  [58] DT_0.28                htmlwidgets_1.6.2      httr_1.4.6            
##  [61] gplots_3.1.3           RColorBrewer_1.1-3     ellipsis_0.3.2        
##  [64] farver_2.1.1           pkgconfig_2.0.3        XML_3.99-0.14         
##  [67] sass_0.4.6             utf8_1.2.3             STRINGdb_2.10.1       
##  [70] labeling_0.4.2         tidyselect_1.2.0       rlang_1.1.1           
##  [73] later_1.3.1            munsell_0.5.0          tools_4.2.3           
##  [76] cachem_1.0.8           cli_3.6.1              gsubfn_0.7            
##  [79] generics_0.1.3         RSQLite_2.3.1          fdrtool_1.2.17        
##  [82] evaluate_0.21          fastmap_1.1.1          mzID_1.36.0           
##  [85] yaml_2.3.7             knitr_1.43             bit64_4.0.5           
##  [88] caTools_1.18.2         KEGGREST_1.38.0        ncdf4_1.21            
##  [91] mime_0.12              compiler_4.2.3         rstudioapi_0.14       
##  [94] plotly_4.10.2          png_0.1-8              affyio_1.68.0         
##  [97] stringi_1.7.12         bslib_0.5.0            highr_0.10            
## [100] MSnbase_2.24.2         lattice_0.21-8         ProtGenerics_1.30.0   
## [103] Matrix_1.5-4.1         tmvtnorm_1.5           vctrs_0.6.3           
## [106] pillar_1.9.0           norm_1.0-11.0          lifecycle_1.0.3       
## [109] BiocManager_1.30.21    jquerylib_0.1.4        MALDIquant_1.22.1     
## [112] GlobalOptions_0.1.2    data.table_1.14.8      cowplot_1.1.1         
## [115] bitops_1.0-7           httpuv_1.6.11          R6_2.5.1              
## [118] pcaMethods_1.90.0      affy_1.76.0            promises_1.2.0.1      
## [121] KernSmooth_2.23-21     codetools_0.2-19       MASS_7.3-60           
## [124] gtools_3.9.4           assertthat_0.2.1       chron_2.3-61          
## [127] proto_1.0.0            rjson_0.2.21           withr_2.5.0           
## [130] GenomeInfoDbData_1.2.9 parallel_4.2.3         hms_1.1.3             
## [133] grid_4.2.3             rmarkdown_2.22         shiny_1.7.4
knitr::knit_exit()